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A method  is  presented  for  modeling  a transmission 
line  for  transient  analysis  study  on  the  digital  computer. 
The  line  is  modeled  as  a finite  number  of  lumped  parameter 
sections.  Each  section  is  modeled  in  an  equivalent  section 
of  resistors  and  current  sources  developed  from  solving  the 
voltage  and  current  equations  by  the  trapezoidal  rule  for 
integration.  The  integration  takes  place  over  a period  of 
time  from  a known  state,  t,  to  an  unknown  state,  t+At.  The 
time  step.  At,  is  taken  to  be  the  lossless  travel  time  for 
the  traveling  wave  to  cross  each  section.  The  single  phase 
lossless  case  is  handled  first,  then  losses  are  accounted 
for,  and  finally  the  three  phase  line  is  dealt  with. 
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I. 


INTRODUCTION 


The  need  to  study  the  transient  phenomena  on  power 
transmission  lines  results  from  the  high  voltages  experi- 
enced on  the  line  due  to  normal  but  abrupt  switching  actions. 
These  high  voltages  that  appear  are  on  an  order  of  magnitude 
of  two  to  three  times  the  rated  line  voltage.  Transients 
can  be  caused  by  other  factors,  such  as  atmospheric  dis- 
turbance, but  the  majority  is  due  to  normal  switching  opera- 
tions on  the  line.  These  transients  usually  last  for  only 
a few  milliseconds  [1] , but  insulators  and  other  equipment 
can  be  permanently  damaged.  The  study  of  transients  on 
transmission  lines  has  been  underway  for  many  years  with 
the  classical  equations  being  well  known.  With  the  advent 
of  the  digital  computer,  methods  are  now  available  to  solve 
the  classical  equations  numerically  and  with  a high  degree 
of  accuracy.  This  thesis  deals  with  modeling  the  trans- 
mission line  for  the  digital  computer  in  order  to  solve  for 
the  transient  voltages  and  currents  that  exist  due  to  switch- 
ing actions  that  can  occur  from  energizing  or  deenergizing 
the  line. 

The  solution  to  the  classical  transmission  line 
equations  is  well  known  [1,  2,  3]  and  need  not  be  presented 
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here.  The  nature  of  the  solutions  results  in  the  concept 
of  traveling  waves  on  the  transmission  line.  Since  the 
transmission  line  parameters  of  resistance,  inductance, 
capacitance,  and  conductance  are  distributed  uniformly 
throughout  the  line,  this  provides  the  line  with  its  wave 
carrying  capability.  It  is  much  like  any  other  physical 
continua,  such  as  air  and  water,  in  this  respect  [1] . 

These  traveling  waves  on  the  line  are  of  two  types — 
forward  traveling  and  reverse  traveling.  The  reverse 
traveling  wave  is  a scaled  version  of  the  forward  travel- 
ing wave.  This  scaling  factor  is  called  the  reflection 
coefficient.  Solutions  have  been  very  complicated  except 
for  the  simplest  cases  and  have  typically  dealt  with  a 
lossless  line,  i.e.,  resistance  and  conductance  are  assumed 
to  be  zero.  One  such  solution  utilizes  the  Bewley  Lattice 
diagram  which  requires  that  the  reflection  coefficients  for 
the  sending  and  receiving  ends  be  calculated  [4]. 

Most  all  the  work  done  in  these  solutions  is  for  a 
single  phase  line.  When  three  phase  lines  are  studied, 
the  concept  of  three  phase  is  lost  because  of  the  transient 
phenomena.  The  line  can  be  viewed  as  three  separate 
phases  by  using  a matrix  transformation  to  decouple  the 
phases.  This  approach  is  used  in  this  thesis.  In  modeling 
the  three  phase  line  the  earth  return  for  the  ground  cur- 
rents must  be  included,  and  in  transient  analysis  this  intro- 
duces the  complex  situation  of  handling  the  frequency  dependency 
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of  resistance  and  inductance  of  the  ground  mode  [5] . This 
topic  will  be  discussed  later. 

It  has  been  mentioned  that  the  transmission  line  is 
composed  of  uniformly  distributed  parameters.  This  thesis 
models  the  line  as  a finite  number  of  sections  each  having 
lumped  parameters,  as  shown  in  figure  1-1.  The  argu- 
ment for  this  is  that  as  the  number  of  sections 
approaches  infinity  as  their  lengths  become  smaller  and 
smaller,  it  approximates  the  distributed  line.  Initially, 
the  line  will  be  considered  lossless  with  the  lossy  case 
being  handled  later. 

In  researching  the  literature  a paper  by  Hermann 

W.  Dommell,  "Digital  Computer  Solution  of  Electromagnetic 

Transients  in  Single-  and  Multiphase  Networks"  [5] , was 
i- 

listed  as  a source  by  nearly  everyone  who  was  working  on 
the  problem  of  digital  solutions  to  transmission  line  tran- 
sients. His  method  for  solving  transients  is  to  handle  the 
distributed  parameters  with  a method  called  characteristics 
\ and  the  lumped  parameters  with  the  trapezoidal  rule  for 

integration.  The  method  of  characteristics  is  described 
more  fully  in  a paper  by  F.  H.  Branin,  Jr.,  "Transient 
Analysis  of  Lossless  Transmission  Lines"  [6] . The  in- 

K 


elusion  of  frequency  dependent  parameters  in  the  problem  is 
presented  by  Alan  Budner  in  his  paper,  "Introduction  of 
Frequency-Dependent  Line  Parameters  into  an  Electromagnetic 


Figure  1-1.  Transmission  line  represented  by  n-sections  and  n+1  nodes 
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Transient  Program"  [7] . This  same  problem  is  also  discussed 
by  J.  K.  Snelson  in  "Propagation  of  Traveling  Waves  on 
Transmission  Lines — Frequency  Dependent  Parameters"  [8] . 

S.  C.  Tripathy  and  N.  D.  Roa  present  a method  in  "A-Stable 
Numerical  Integration  Method  for  Transmission  System  Tran- 
sients" [9]  for  handling  nonlinear  elements  with  a non- 
iterative technique. 

The  basis  for  this  thesis  is  Dommell's  work.  A method 
will  be  developed  to  handle  the  transmission  line  on  a 
digital  computer  for  transient  analysis.  Lumped  parameters 
for  the  transmission  line  will  be  dealt  with  exclusively. 

The  lossless  single-phase  line  will  be  developed  first,  then 
losses  will  be  accounted  for,  and  finally  the  three  phase 
line  will  be  analyzed. 


II 


THE  SINGLE-PHASE  LINE 


Lossless  Case 

Figure  2-1  presents  a typical  section  of  the  trans- 
mission line  presented  in  figure  1-1.  As  the  line  is 
divided  into  sections,  it  will  have  n sections  and 
n + 1 nodes.  The  development  that  follows  will  lend 
itself  to  digital  computer  techniques.  Since  the  digital 
computer  cannot  give  the  entire  listing  of  a transient  on  a 
transmission  line  [5] , the  development  will  be  one  that 
recognizes  that  it  can  give  the  results  of  computations  at 
some  time  t + At  where  the  results  at  time  t are  known. 
Referring  to  figure  2-1,  the  equation  for  the  current  through 
the  inductor  can  be  written  as, 


di . 

V.  - V.  , , = L — 

i l+l  dt 


dii  " L <Vi  - Vi+1>  dt 


(2-la) 

(2-lb) 


Integrating  from  the  known  state,  t,  to  an  unknown  state, 
t + At,  using  the  trapezoidal  rule  for  integration  [10] , gives 


/•t+At 
! dii 


c t+At 


Vi+1>  dt 


(2-lc) 


6 


8 


, < t 


/ 


v. 


i t -:~ 


At  * J 

ii(t+At)  - ii ( t)  = 2l  lvi(t+At)  - vi+1(t+At) 


f r 


+ v.(t)  - v.+1(t)] 


(2-ld) 


At 


ii(t+At)  = 2l  tvi(t+At)  - vi+1(t+At)] 


At 


+ 2L  [vi(t)  ' Vi+l(t)]  + ii(t) 


(2-le) 


In  equation  2-le  the  current  at  the  ith  node  at  t+At  is 
dependent  on  the  difference  in  voltage  at  the  i and  i+1  nodes 


T 

divided  by  an  equivalent  resistance  between  the  two  nodes, 


The  voltage  and  current  in  equation  2-le  at  time  t can  be 
viewed  as  the  past  voltage  and  current,  and  therefore  are 
known.  Let, 


Ii(t) 


JE  [vi<t)  - vi+l(t)I  + ii(t) 


( 2-lf ) 


These  known  values  at  time  t will  be  viev/ed  as  a current 


source,  1^ ( t) . 


Turning  to  the  capacitor  in  the  section,  the  current 
through  it  is  given  by, 

dv . 


i = C A*1 

cl  L at 


(2-2a) 


The  current  can  also  be  expressed  as. 


1ci  xi  " xi+l 


(2-2b) 


j 
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Substituting  equation  2-2b  into  2-2a  and  rewriting. 


dvi-n  = h !ii  - ii+i)dt 


(2-2c) 


Using  the  trapezoidal  rule  and  integrating  from  t to  t+At, 


t+At 

dv 


i+1 


1 

C 


t+At 


(ii  - ii+!)dt 


(2-2d) 


v.+1(t+At)  - vi+1(t)  = [ij.  (t+At)  - ii+1(t+At) 


+ ii(t)  - ii+1(t)]  (2-2e) 


vi+l (t+At)  = [ii(t+At)  - ii+1(t+At)] 


+ 2C  (t)  ~ ii+i (t)  1 


vi+i(t) 


( 2 — 2 f ) 


The  known  values  in  equation  2-2f  now  appear  as  a voltage 
source.  Let, 


Vi+l(t) 


2C 


ii+l(t)] 


Vi+l(t) 


(2-2g) 


The  equivalent  circuit  for  the  line  section  is  shown  in 
figure  2-2.  Since  the  network  now  lends  itself  to  general 
nodal  analysis, it  is  desirable  to  transform  the  voltage 
source  to  a current  source.  Using  Norton's  Theorem  to 
accomplish  this,  the  resulting  circuit  is  shown  in  figure  2-3 . 


Figure  2-2.  Typical  lossless  transmission  line  equivalent 
section 


2L 

At 


ii+l ( t+At) 


Figure  2 


-3.  Typical  transmission  line  equivalent  section 
with  all  current  sources 
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Starting  with  the  typical  section  of  the  transmission  line 
composed  of  inductance  and  capacitance,  the  equivalent  cir- 
cuit is  one  composed  of  resistive  elements  and  current 
sources.  This  allows  analysis  of  each  section  to  be  done  by 
nodal  techniques  and  without  solving  differential  equations. 
The  entire  equivalent  transmission  line  is  shown  in  figure 
2-4.  For  nodal  analysis , it  is  more  convenient  to  deal  with 
conductance  than  resistance.  In  figure  2-4. 


Gs 

At 

2L 

and, 

(2-3) 

G 

P 

2C 

At 

(2-4) 

Writing  a matrix  equation  for  the  entire  line  using 
conventional  nodal  analysis, 

[Y]v  = C ( 2-5a) 


where,  excluding  the  end  nodes, 


2G  + G , i=j=2,3,  ...,n 

s p 

-Gs  , i = j+1  ( 2-5b) 

0 , i = j+2 , 3 , . . . , n 


The  general  entry  for  v is  v^(t+At)  and  for  C,  a current 
vector,  is  given  by, 


Figure  2-4.  Lossless  transmission  line 
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ci  = Ii_1<t)  + Gp’Vi<t>  i = 2,  . . . , n (2-5c) 

where,  again  excluding  the  end  nodes, 


Ii-l(t)  = Gs[vi-l(t)  ' vi(t)l  + 


Vi(t)  * + Vi(t) 


- Gs[vi(t)  - vi+1(t)  ] i (t) 


The  solution  to  equation  2-5a  is  given  by, 


v = [Y]  1C 


( 2-5d) 


(2-5e) 


(2-5f ) 


(2-5g) 


A method  for  inverting  the  Y-matrix,  which  is  a sparse 
matrix,  is  given  in  Appendix  A. 

The  source  for  the  line  was  chosen  to  be  a current 
source  with  shunted  inductance,  capacitance,  and  resistance 
as  shown  in  figure  2-5.  Each  element  of  the  source  will  be 
handled  separately  in  order  to  develop  a model  compatible 
with  the  line  model. 

For  the  inductor. 


di. 


v.  = L 
1 s dt 


( 2-6a) 


Using  the  trapezoidal  rule  and  integrating  from  t to  t+At, 


ft+At 


ft+At 


diL  " !T 
s s 


v^dt 


(2-6b) 
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iLs<t+At)  - i (t)  = [vx(t+At)  + v1  <t)  ] ( 2-6c ) 

iL  (t+At)  = 4n~  Vl(t+At>  + V;L(t)  +iL  (t)  (2-6d) 


!l  (t)  = 


2IT  Vfc>  ~ iL  (t) 

s s 


(2-6e) 


then. 


iL  (t+At>  = vx (t+At)  - I (t)  ( 2-6 f ) 


Similarly,  for  the  capacitor. 


iCs  Cs  dt 


2Cs  ^ 


ic  (t+At)  = V]L(t+At)  - V;L(t)  - 


(2-7a) 


ic  (t)  (2-7b) 

s 


IC  (t)  = "ST  Vl(t)  + iC  (t)  ( 2-7c) 


then. 


iC  (t+At>  = vx(t+At)  - Ic  (t)  (2-1  d) 

s s 

Since  the  resistor  does  not  contribute  a current  source 
in  modeling,  it  remains  unchanged.  The  equivalent  source 
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model  is  shown  in  figure  2-6.  In  order  to  simplify  the  model 
and  make  it  more  compatible  with  the  line  and  also  for 
nodal  analysis,  the  current  sources  are  combined  and 
resistive  elements  are  also  combined  but  as  conductances. 
The  final  result  is  shown  in  figure  2-7,  with, 


y<=>  = y«  + !l  <«  + iL 

s s 


(2-8) 


r - At  s , 1 

g 2L  At  + R 
3 s s 


(2-9) 


The  receiving  end  termination  can  also  be  modeled 
in  a general  circuit  of  shunt  resistance,  inductance  and 
capacitance.  Following  the  same  argument  as  before  in  the 
source,  the  receiving  end  circuit  for  a generalized  load  is 
shown  in  figure  2-8.  As  before. 


G - At  2CL  1 
GL  2L  At  Rt 

J-i  ij 


(2-10) 


IL(t)  = Ic  (t)  + IL  (t) 
L L 


( 2-lla) 


where. 


IC  (t)  At  vn+l (t)  + XC  (t) 
J_»  L 


(2-llb) 


^■L.  (t)  2L  Vn+l(t)  1L  (t) 

Jl»  L 


(2-llc) 


Returning  to  the  nodal  equation. 


Figure  2-6.  Equivalent  source  model 


g 
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[Y]v 

= c 

(2-5a) 

the  entries 

of  the  Y-matrix  and 

C-vector  can  now  be 

com- 

pleted.  For 

the  Y-matrix, 

yn  ' G, 

+ G 

3 s 

(2-5h) 

where. 

r _ At 

Gg  2L 

2C 

h s + 1 

At  R 

s 

(2-9) 

G 

s 

At 

2L 

(2-3) 

and. 

yn+l , n+1  - Gs 

+ G + G. 

P L 

(2-5i) 

where , 

GP  ’ 

2C 

At 

(2-4) 

For  the  C-vector, 

°i = Vfc) 

- I1(t) 

( 2-5  j ) 

where , 

Vfc)  = Is(t)  + I 

: (t)  + i (t) 

S S 

(2-8) 

Ix(t) 

= (t)  - v2  (t)  ] 

+ i-^t) 

(2-5f ) 

and. 

cn+l 

' V**  * GpVn+l(t) 

+ lL(t) 

(2-5k) 

where. 
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Xn(t)  = Gs[vn(t)  " vn+l(t)]  + in(t)  (2"5f) 


G V - (t)  = i (t)  - i . (t)  + G v . (t)  (2-5e) 

p n+1  n n+1  p n+1 


IL(t) 


lr  (t) 


+ IT  (t) 


(2-lla) 


Lossless  Case  Examples 

An  example  problem  was  chosen  to  apply  to  the  pre- 
ceding development.  Computer  input  data  for  the  cases  tested, 
2,  10,  and  40  section  lines,  are  shown  in  figure  2-9,  with 
only  the  number  of  sections  changing  for  each  example. 

Data  were  chosen  to  facilitate  manual  calculations.  The 
velocity  of  wave  propagation  on  a lossless  line  is  given  by 
[11]  , 

velocity  = — (2-12) 

/LC 


and  the  characteristic  impedance  is  given  by  [11] , 


Z = 
o 


(2-13) 


For  the  example  chosen,  both  of  these  are  one  and  the  sending 
and  receiving  ends  are  terminated  in  the  characteristic 
impedance.  With  this  configuration,  there  should  be  no 
reflected  voltages  along  the  line  and  with  the  one  amp  cur- 
rent source,  voltage  and  current  should  stabilize  at  the 
same  value,  0.5.  The  delta  t chosen  for  the  equations 
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Figure  2-9.  Computer  input  for  2,  10,  and  40  section  line  examples 
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corresponds  to  the  travel  time  for  the  line  section.  Results 
for  the  three  cases  tested,  2,  10,  and  40  sections,  are 
shown  in  figures  2-10,  2-11,  and  2-12,  respectively.  The  plots 
show  voltage  as  function  of  position  at  a time  that  corres- 
ponds to  the  wave  traveling  halfway  down  the  line.  The 
theoretical  wave  shapes  are  shown  as  dashed  lines.  As  might 
be  expected,  the  40  section  line  exhibited  a waveform  that 
more  closely  approximated  the  theoretical,  and  is  more 
oscillatory  in  nature  than  the  other  cases.  The  calcula- 
tions for  the  three  cases  were  allowed  to  continue  for  a 
time  period  until  they  stabilized  and  these  results  are 
shown  in  figures  2-13,  2-14,  and  2-15. 

In  order  to  see  the  effect  of  the  size  of  delta  t, 
a smaller  delta  t than  for  any  previous  case  was  chosen, 

At  = .01,  and  the  three  cases  run  again.  These  results  are 
shown  in  figures  2-16,  2-17,  and  2-18  for  a travel  time  of  halfway 
down  the  line.  When  compared  with  figures  2-10  , 2-11 , and  2-12 , 
respectively,  the  wave  shapes  do  not  appear  very  different, 
except  that  they  are  generally  steeper  at  the  leading  edge. 

To  further  check  its  effect,  runs  were  made  with  the  40  sec- 
tion line  for  varying  delta  t’s,  larger  and  smaller  than  the 
travel  time  for  each  section.  These  results  are  shown  in 
figures  2-19,  2-20,  and  2-21.  The  smallest  delta  t chosen  was 
0.001,  figure  2-21,  and  results  do  not  vary  appreciably  from  that 
of  0.005  in  figure  2-20.  These  results  confirm  Dommell's 
results  [5]  that  changing  delta  t tends  to  change  the  phase 


Figure  2-10.  Voltage  versus  position  for  a 2 section  line 


Figure  2-11.  Voltage  versus  position  for  a 10  section  line,  At 


-14.  Voltage  versus  position  for  a 10  section  line.  At 


Figure  2-15.  Voltage  versus  position  for  a 40  section,  At 


Figure  2-17.  Voltage  versus  position  for  a 10  section  line.  At  - .01 
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sition  for  a 40  section  line 


position  for  a 40  section  line,  At 


Figure  2-21.  Voltage  versus  position  for  a 40  section  line.  At  = .001 
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position  of  the  high  frequency  oscillations  but  not  their 
amplitudes.  This  does  not  make  the  choice  of  delta  t cri- 
tical and  for  the  remainder  of  this  thesis  delta  t will  be 
chosen  to  be  the  lossless  travel  time  for  the  line  section. 
The  length  of  the  line  section  was  determined  from  the  fact 
that  most  delta  t's  are  in  the  neighborhood  of  50  micro- 
seconds. Since  the  velocity  of  a wave  on  a lossless  line  is 
approximately  the  speed  of  light,  the  wave  travels  approxi- 
mately 15  kilometers  (9.3  miles)  in  50  microseconds.  This 
figure  of  15  kilometers  is  used  to  determine,  to  the  nearest 
whole  number,  the  number  of  sections  needed  to  represent 
the  line  under  consideration. 


! 
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t+At 


di.  = f 
l L 


t+At 


(Vi  " Rii  ' Vi+l)dt 


(2-14c) 


ii  (t+At)  - i.(t>  = §i{v.(t+At)  + v.(t)  - Rfi.(t+At) 


+ ii(t) ] - Ivi+1(t+At) 


+ vi+1(t)]} 


[i  + (t+At)  - |i[Vi(t+At)  - V (t+At)] 


i+1 


At 


+ 2L  [vi(t)  " vi+i(t)] 


, r,  RAt,  . 

+ [i  - 


1i(t+At)  - t 2L  1 ] [v±(t+At)  - v.  (t+At)] 

~ + R 1 1+1 


At 


+ [ 


2L—  ] [v.(t)  - vi+1(tn 


At 


2L  _ 

At  R 


2L 


At 


+ R 


ii(t) 


Let, 


Vo 


= 1 STT  1 lvi(tl  " W*” 


At 


2L 

At  R 


2L 


At 


+ R 


( t) 


(2-14d) 


( 2-14e) 


(2-14  f) 


(2-14g) 
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The  model  for  the  remainder  of  the  section  remains  the  same 
as  the  lossless  model.  The  lossy  circuit  model  appears 
now  as  in  figure  2-23,  where, 

vi+1(t)  - §§  [ii(t)  - ii+1(t)l  + vi+1(t>  (2-2g) 

If  resistance  is  set  to  zero  in  this  model,  then 
it  reduces  to  the  lossless  case.  The  same  nodal 
equations  can  be  written  as  in  the  lossless  case  with  the 
appropriate  changes  made  to  the  Y-matrix  and  C-vector. 

Lossy  Example 

An  example  was  chosen  from  Dommell's  work  [5]  to 
compare  results  of  the  two  programs.  The  line  data  is  listed 
below. 

line  length  = 320  miles 

R = 0.0376  ohms/mi 

L = 1.52  mH/mi 

C = 0.0143  pF/mi 

line  termination  = 0.1  H 

source  = 10  V.  step 

The  data  was  scaled  in  order  to  input  it  to  the  computer. 

The  scaling  method  is  described  in  Appendix  B.  Since  Dom- 
mell's program  used  32  sections  to  represent  the  line,  the 


same  number  was  used  in  this  example.  The  results  are  shown 
in  figure  2-24.  The  agreement  between  the  two  program  results 


1 i l t+AtJ 


i-j+1  (t+At) 


-24 


4 3 


is  very  good.  Wave  shape  is  almost  identical  with  slight 
variation  in  amplitudes.  This  possibly  is  due  to  the  fact 
that  Dommell's  method  for  handling  resistance  is  different 
from  this  thesis'  method.  Also,  it  is  not  precisely  known 
how  Dommell  handled  end  effects. 

Modeling  a Sinusoidal 
Voltage  Source 

Modeling  an  ideal,  non-time  varying  voltage  source 
as  a current  source,  as  in  the  previous  example,  proved  to 
be  straight  forward.  Technically,  the  Norton  Equiva- 
lent of  an  ideal  voltage  source  is  an  infinite  current 
source  shunted  by  a zero  resistance.  For  the  computer  pro- 
gram developed  for  this  analysis,  the  ideal  voltage  source 
was  modeled  as  a very  large  current  source  shunted  by  a 
very  small  resistance.  The  values  were  chosen  such  that  the 
open  circuit  voltage  equaled  that  of  the  voltage  source. 

This  method  also  works  for  a sinusoidal  source  but  it  is 
modified  slightly. 

When  a sinusoidal  source  is  used  as  a model  for 
energizing  a transmission  line  it  is  usually  modeled  as  a 
generator  with  a series  impedance  made  up  of  inductance  and 
and  resistance.  Choosing  a model  for  the  network  that 
energizes  a transmission  line  is  a non-trivial  task,  but  it' 
is  not  the  sub ject  of  this  thesis . In  this  thesis , a sinusoidal 
voltage  source  with  series  impedance  is  modeled  in  the  following 
way.  The  source  itself  is  always  modeled  as  an  ideal 
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voltage  source,  while  the  series  impedance  is  added  as  a 
new  section  to  the  beginning  of  the  line.  This  new  sec- 
tion is  handled  as  was  the  lossy  section  handled  earlier  in 
this  chapter  (Equations  2-14a-g) . This  method  increases 
the  number  of  sections  and  nodes  by  one.  This  development 
becomes  more  important  when  three  phase  lines  are  encoun- 
tered. The  equivalent  model  for  a general  sinusoidal  volt- 
age is  shown  in  figure  2-25. 

When  this  source  model  is  used  to  energize  the  line, 
the  Y-matrix  and  C-vector  are  modified  as  follow.  Refer- 
ring to  figure  2-25,  let 


G = 1/[(2L  /At)  + R ] 
gs  s s 


(2-15) 


Equation  2-9  now  becomes 


G = 10 
9 


(2-9'  ) 


Then 


yn  0 - G + G 
0 , 0 g gc 


( 2-5it ) 


and , 


y0,l  - yl,0  - -Gg. 


(2-5ra; 


Equation  2-5h  now  becomes 


Yn  = G + G 

11  s 

s 


(2-5h ' ) 
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and  equation  2-5b  becomes 


2G  + G 
s p 


0 


* i = j / j 
, i = j+1,  j 
9 i = j^2 t 3 / • 


2/3/  . • • # 
2/3/  • • . / 

• • z n 


except  as  noted.  Equation  2-5i  remains  unchanged 
C-vector  changes,  let 

IQ(t)  = [ 2l  ][vQ(t)  - v1(t)] 

s . r. 


i0(fc) 


Equation  2-8  now  becomes 


2L 
s 

At 


- R 


2L 

+ R 

At  S 


I (t)  = 106 [E  Cos (wt  + 0)  ] 
g m 


then , 


c 0 - 

Equation  2-5j  now  becomes 

c1  = lQ(t)  - I1(t) 

and  equation  2-5c  is  valid  except  the  subscript  i 


( 2-5b 1 ) 


For  the 


( 2-5n) 


(2-8'  ) 


( 2-5p) 


( 2 - 5 j ' ) 

now  starts 


at  3 instead  of  2. 


III.  THE  THREE  PHASE  LINE 


k 


Three  Phase  Model 

In  dealing  with  the  three  phase  line  under  transient 
conditions,  it  is  desirable  to  analyze  each  phase  separately. 
To  accomplish  this,  the  phases  must  be  decoupled  because  of 
the  mutual  inductances  that  exist  between  them.  This  can  be 
done  with  a similarity  transformation  matrix,  known  as  a 
modal  transformation  matrix,  that  diagonalizes  the  line 
impedance  and  admittance  matrices  12].  In  Dommell's  paper 
[6] , he  introduces  a modal  transformation  matrix  [T]  for  a 
three  phase  line  as, 


with 


1 


[T]-1  = 1/3 


1 


1 

-1 


1 

0 


10-1 


(3-1) 


This  matrix  is  only  valid  for  a completely  transposed  line. 
There  are  other  modal  matrices  which  are  used  on  all  types 

of  lines  but  they  will  produce  off-diagonal  elements  in  the 
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transformed  matrices.  These  off-diagonal  elements  are 

1 

small  when  compared  with  the  diagonal  elements,  and  are 
generally  ignored.  This  is  strictly  true  assuming  a 
totally  transposed  line,  which  will  not  produce  off- 
diagonal  elements  when  transformed.  In  this  work,  off- 
diagonal  elements  will  be  ignored.  This  allows  the  sequence 
values  to  be  substituted  for  the  modal  values  since  they 
are  equal.  Specifically 
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The  phase  voltages  and  currents  are  defined  in 
terms  of  the  modal  values  as  follows: 


abc 


= [T]v 


oag 


(3-3) 


and. 


[T]  i 


oag 


(3-4) 


With  the  line  defined  now  with  its  modal  values,  the 
problem  reverts  back  to  the  single  phase  case  as  des- 
cribed in  Chapter  2.  Each  mode  will  be  treated  as  the 
equivalent  line  in  figure  2-4,  except  with  losses.  After 
each  mode  is  solved,  essentially  three  single  phase 
problems,  equations  3-3  and  3-4  will  be  used  to  find 
the  phase  values. 

The  end  effects  for  the  three  phase  line  are 
essentially  handled  as  in  the  single  phase  case.  The 
three  phase  network  energizing  the  line  will  be  assumed 
to  be  a three  phase  voltage  source  with  series  imped- 
ance. The  voltage  source  will  be  handled  as  described 


in  Chapter  2. 
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Three  Phase  Example 

A three  phase  example  problem  was  chosen  from  work 
done  by  Southern  Company  Services,  Inc.,  with  their  tran- 
sient program  "Surge."  The  line  and  system  data  are  listed 
below. 


System:  345  KV,  100  MVA,  50  Hz 


Source  data 

positive  sequence  voltage  (p.u.): 


1.001l[_0°  (line  to  neu- 
tral peak  value) 


impedance  (p.u.):  0.0115  + j0.2206 

switching  angles:  A = 71.8°  (3.99  ms) 

B = 163.1°  (9.06  ms) 

C = 32.0°  (1.78  ms) 

(note:  Switching  angles  are  used  to  simulate  assyn- 

chronous  switching.) 


Line  data 


Load  data 


zero  sequence: 


pos/neg  sequence: 


length: 


R = 0.418  ohm/mi 
L = 5.198  mh/mi 
C = 0.01232  yf/mi 

R = 0.0644  ohm/mi 
L = 1.629  mh/mi 
C = 0.01908  yf/mi 

126  mi 


open  circuit 


The  results  from  this  example  are  shown  in  the  com- 
puter plot  of  the  receiving  end  voltage  in  figure  3-1.  The 
waveforms  and  amplitudes  are  in  excellent  agreement  with 
"Surge"  results.  The  "Surge"  results  have  slightly  lower 
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voltage  maximum  but  that  is  to  be  expected  since  it  handles 
the  frequency-dependence  of  parameters  which  is  discussed 
in  the  next  section. 

Frequency-Dependence  of 
Line  Parameters 

This  thesis  has  not  attempted  to  include  a method  for 
handling  frequency-dependent  line  parameters.  Other  works 
in  this  area  have  dealt  with  the  problem,  and  since  it 
does  have  a bearing  on  the  transients  observed  on  the  line 
it  will  be  mentioned  here. 

An  overhead  transmission  line  is  composed  of  a cer- 
tain number  of  phase  conductors  and  neutrals.  The  phase 
conductors  are  separated  from  each  other  but  the  neutrals 
are  connected  through  the  towers  and  are  thus  grounded.  The 
modeling  of  this  ground  return  for  inclusion  into  a 
model  of  the  line  is  very  complex  due  to  the  non- 
uniformity of  the  earth.  In  an  early  paper  by  Carson  on 
this  subject  [14],  he  established  the  fact  that  for  a single 
conductor  with  ground  return  its  resistance  and  inductance 
per  unit  length  are  proportional  to  frequency  (f  in  Hz)  . Another 
author  [8]  has  noted  Carson's  results  as: 

R « (f)k  (3-5a) 


where. 


0.5  < k < 1.0 


( 3-5b) 


L « (f)n  (3-6a) 


where. 


-0.5  < n < 0 ( 3-6b) 


In  three  phase  transient  analysis  using  modal  tech- 
niques, the  O-mode  is  often  referred  to  as  the  ground  mode. 
In  a study  of  frequency  effects  on  modal  values  by  Hedman 
[5]  , it  was  found  that  the  mode  most  affected  by  frequency  is 
the  ground  mode.  His  conclusions  on  earth  affects  are 
listed  below. 

1.  Carson's  earth-correction  terms  produce  the  predomi- 
nant earth-correction  effects  for  a transmission 
line  over  an  imperfect  earth. 

2.  Carson's  earth-correction  resistance  terms  are 
proportional  to  frequency  and  to  the  square  root 
of  frequency,  respectively  in  the  low-  and  high- 
frequency  regions. 

3.  Effects  of  the  high  relative-dielectric  constant  of 
the  earth  are  significant  only  for  frequencies 
higher  than  0.5  MHz  and  when  both  earth  resistivity 
and  dielectric  constant  are  high. 

4.  Earth  correction  for  admittance  terms  appear  to  be 
unimportant  for  frequencies  lower  than  1 MHz. 

5.  Carson's  earth  correction  terms  significantly  affect 
the  modal  voltages  and  eigenvectors  for  frequencies 
from  60  Hz  to  1 MHz. 

6.  Modal  analysis,  using  the  perfect  earth,  should  be 
adequate  for  radio-noise  propagation  studies. 


7.  For  carrier-current  analysis,  earth  effects  become 
significant. 
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The  frequency  dependent  resistance  and  inductance  have  a 
damping  effect  on  transient  voltages  when  compared  to  tran- 
sients that  do  not  consider  frequency  dependence. 

In  computer  programs  that  deal  with  frequency 
dependence  a frequency  domain  technique  is  used  to  deter- 
mine the  values  of  parameters  in  the  equations  already 
present.  Methods  such  as  the  Fourier  Transform  [8,  15]  and 
the  Modified  Fourier  Transform  [9],  are  used  to  evaluate  the 
parameters  over  a range  of  frequencies  at  each  time  step  in 
the  program.  A typical  range  of  frequencies  would  be 
0-12.8  kHz  [8]  . 


IV.  CONCLUSION 

The  method  for  modeling  the  transmission  line  that 
results  from  the  trapezoidal  rule  of  integration  is  a very 
straight  forward  way  of  solving  transmission  line  transients. 
In  fact,  at  each  time  step  the  problem  to  be  solved  is  that 
of  a d.c.  circuits  problem.  The  argument  for  this  is  that 
the  time  step  is  selected  small  enough,  i.e.  the  lossless 
travel  time  for  the  traveling  wave  to  cross  the  section, 
that  nothing  changes  during  that  span  of  time.  Handling 
the  end  effects  of  the  line  using  the  trapezoidal  rule 
proved  to  be  very  compatible  with  the  rest  of  the  line 
model . 

The  argument  was  made  that  as  the  number  of  sections 
increases  while  each  section  length  decreases  that  this  more 
closely  approximated  the  actual  line  performance.  This 
proved  to  be  true  in  the  sample  lines  of  2,  10,  and  40  sec- 
tions in  the  single  phase  case.  The  number  of  sections  and 
the  time  step  chosen  were  related  by  selecting  a time  step 
that  would  equal  the  lossless  travel  time  for  the  section. 

A smaller  time  step  proved  to  be  mere  accurate  for  the 
smaller  number  of  sectioned  lines  but  about  the  same  for 
the  larger  number  of  sectioned  lines. 


I 
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The  modal  technique  used  in  the  three  phase  case 
proved  to  be  very  powerful  in  handling  three  phase  tran- 
sients. Its  decoupling  of  the  phases  into  the  modal  values 
just  presented  the  problem  of  solving  three  single  phase  i 

cases.  Transformation  back  to  the  phase  values  presented  1 

the  desired  results.  Although  frequency-dependence  of  line 
parameters  was  not  included,  it  did  not  present  a serious 
problem  in  the  analysis.  The  analysis,  as  developed,  pro-  : 

duces  slightly  higher  voltages  than  had  frequency-dependence 
been  included.  Since  the  maximum  voltages  are  of  prime  j 

interest  in  transient  analysis,  this  places  this  thesis' 
results  on  the  conservative  side  in  determining  them. 

Finally,  the  examples  cited  and  run  on  the  program 
developed  present  excellent  agreement  between  this  method 
and  the  methods  previously  developed. 
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INVERSION  TECHNIQUE  FOR  A SPARSE  MATRIX 

The  general  form  of  the  Y-matrix  discussed  in  Chap- 
ter 2 is  shown  in  figure  A-l.  In  a typical  computer  routine 
for  inverting  a matrix,  every  entry  would  be  used  in  deter- 
mining the  inverse.  Since  the  Y-matrix  is  a tridiagonal 
matrix,  it  would  be  advantageous  to  exploit  its  sparseness 
for  the  computer. 

The  method  chosen  for  inverting  the  matrix  is  the 
Gauss-Jordon  method  [12]  . This  method  uses  an  augmented 
matrix  composed  of  the  matrix  to  be  inverted  and  an  identi- 
cal sized  identity  matrix,  as  illustrated  below. 

[A  J I]  (A-la) 

Row  and/or  column  operations  are  performed  on  the  matrix 
to  be  inverted  while  the  same  operations  are  performed  on 
the  identity  matrix.  When  the  A-matrix  has  been  reduced 
to  the  identity  matrix,  the  right  side  of  the  augmented 
matrix  now  contains  A ^ as  shown  below. 

[I  ! A-x]  ( A-lb) 

Since  the  typical  row  entry  of  the  Y-matrix  only 

contains  elements  in  the  y. . and  y.  . 

•'ll  J i , i±l 


positions,  see 


1 
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figure  A-l,  only  these  positions  are  dealt  with  in  the  com- 
puter routine.  Also,  in  the  computer  routine  developed  only 
the  y^  entries  are  changed,  while  all  operations  that 
would  normally  be  performed  on  the  Y-matrix  are  done  only 
on  the  identity  matrix.  The  new  diagonal  entry  is  given  by, 


yi+l,i+l  yi+l,i+l"  (yi+l,i/yi,i)yi,i+l'1-1'  * * * ,n_1 


(A-2 ) 


As  these  operations  are  performed,  the  Y-matrix  is  changed 
to  upper-triangular  form  as  shown  below. 


yll 

y12 

0 

0 

0 . . 

. 0 

0 

y22 

y2  3 

0 

0 . . 

. 0 

0 

0 

y33 

y34 

0 . . 

. 0 

To  further  reduce  Y to  a diagonal  matrix,  row  operations 
are  performed  to  eliminate  the  off  diagonal  entries.  Again 
operations  are  performed  only  on  the  identity  matrix. 
Realizing  that  only  the  diagonal  elements  o ne  Y-matrix 
need  to  be  changed,  since  all  other  elements  are  eliminated, 
saves  the  programmer  and  computer  time.  Now  that  the  Y- 
matrix  is  in  diagonal  form,  each  row  of  the  identity  matrix 
is  divided  by  the  appropriate  yj^  element.  The  inverse  of 
the  Y-matrix  is  now  formed  in  the  place  of  the  identity 
matrix  and  the  Y-matrix  is  set  equal  to  it. 
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SCALING  OF  DATA  FOR  COMPUTER  INPUT 

In  order  to  avoid  working  with  very  small  and  very 
large  numbers  associated  with  a transmission  line,  a 
scaling  method  was  devised  to  input  data  into  the  computer. 
Starting  with  the  transmission  line  equations  (primes  denote 
per  unit  length) , 


+ L' 
a 


(B-la) 


C' 

a 


(B-2a) 


the  bases  for  the  individual  values  are  chosen.  The  sub- 
script a,  denotes  actual  (SI)  values.  Let 


x 

X = - — — (B-3a) 

xbase 


with 


Xbase  = line  len9th  = d (B-3b) 

Let 

L' 

L = v,  9 = 1 (B-4a) 

■^base 


Lbase 


where 


L' 

a 


(B-4b) 


66 


Let 


where 


C' 

C = a - i 
C'  ~ 1 
base 


C'  = r ' 
base  a 


(B-5a) 


(B-5b) 


In  three  phase  analysis,  L£ase  and  C^ase  are  chosen  to  be 
the  positive  sequence  LJ  and  . Time  is  scaled  as 


t = 


'base 


{B-6 a) 


where 


fcbase  ” lossless  travel  time  for  the  line 


t =/L 1 C ' d 

base  a a 


(B-6b) 


Again,  in  three  phase  analysis,  ^ and  are  chosen  tobe 

the  positive  sequence  LJ  and  CJ.  The  voltage  is  scaled  as 


v 


v = 


base 


(B-7a) 


where 


Vbase  " VLN  (line  to  neutral  rated  maximum)  (B-7b) 


To  more  clearly  illustrate  vbase,  assume  a 500  KV  system. 
Then, 


w ^ KV 


Let 


zbase  lossless  characteristic  impedance 


where 
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"base 


C' 

a 


(B-8 ) 


Let 


base 


Vbase^Zbase 


(B-9) 


Substituting  these  values  into  equation  B-la  and  B-2a, 


3 (w,  ) 

base 

3 (xd) 


= (R*  + LL' ) ii, 

a base  a^d/L/C1")  base 

a a 


(B-lb) 


and, 


3 (ii,  ) 
base 

3 (xd) 


= CC 


3 ( vv.  ) 

base 

base  Sftd/lTcV) 
a a 


(B-2b) 


Clearing  terms  on  each  side  of  equations  B-lb  and  B-2b, 


g » * L JL)i 

3X  zbase  5t 


(B-lc) 


and. 


ii  _ P 9v 
3x  ^ 3 1 


(B-2c) 


Letting, 


R = 


R' 

a 


"base 


(B-ld) 


equation  B-lc  becomes 
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The  only  parameter  now  left  to  scale  is  frequency.  Fre- 
quency is  scaled  by  keeping 


wt  = to  t 
a a 


But  time  has  already  been  scaled  and 


to 

a 


to 

base 


t 

t 


a 

base 


cot 


(B-lOa) 


(B-lOb) 


to 


base 


'base 


/L'C'd 
a a 


(B-lOc) 


For  scaling  inductance  (£)  and  capacitance  (c)  that  are  not 
psr  unit  length,  the  following  method  is  used. 


R 

a 


+ jto  Z 

a 

zbase 


a 


R + jto£ 


(B-lla) 


co  z 


a a 


a 


to 

a 


Sl 


1 C ' d 
a a 


Z 


(B-llb) 


* = (L'd)A 

a d. 


(B-llc) 


*’•  *base  = Lad  (B-lld) 

Similarly, 
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An  example  will  help  to  clarify  the  scaling  method.  Using 
the  three  phase  example  of  Chapter  3,  all  values  will  be 
scaled  as  follows: 

Base  data 


xbase 

— 

126  mi 

= d 

L' 

base 

= 

1.629  mh/mi 

C’ 

case 

= 

0.01908 

pf/mi 

^ase 

= 

/LjCjd 

= .7024  ms 

vbase 

= 

345/2 

/3 

= 282  KV 

Ha 

zbase 

= / 

/ i. 

ci 

292  ohms 

CO. 

base 

= 

1//tbase 

= 1423.69  rad/s 

^base 

= 

Lid  = 

0.2053  H 

Cbase 

= 

C^d  = 

2.4041  pf 

Source  data 

Em  = 1.0011 

u)  = [ 2 7T ( 50 ) 1/1423.69  = 0.2207 

0=0° 

R = [0.0115(345) 2/100 ] /292  = 0.0469 
L=  [ (0. 2206 (345) 2/100) 2 (50)1/0.2053=  4.0710 
switching  angles  are  converted  to  times) 

TA  = (71. 8°/360°) (1/50)/. 7024 (10~3)  = 5.6789 


I 


(note : 
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Tb  = (163. l°/360°) (1/50)/. 7024 (10-3)  = 12.9002 
Tc  = (32°/360°) (1/50)/. 7024 (lO-3)  = 2.5310 

Line  data 


zero  sequence:  R = 0. 418 (126) /292  = 0.1804 

L = 5.198(10"3)/1.629(10"3)  = 3.1909 
C = 0. 01232 (10"6)/0. 01908 (10~6)  = 0.6457 

pos/neg  sequence:  R = 0 . 0644 (126 ) /292  = 0.0278 

L = 1.629  (10~3)/1. 629  (10"3)  = 1.0 
C = 0.01908(10"6)/0.01908(10~6)  = 1.0 


Load  data 


open  circuit 


The  bases  for  scaling  are  summarized  in  Table  B-l. 


- 


r 

3 
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VLN(peak) 


V //L' /C' 
VLNX  * a'  a 
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FORTRAN  COMPUTER  PROGRAM 


User's  Guide 

In  order  to  get  data  into  the  single  and  three  phase 
computer  programs,  TTL,  it  must  be  scaled  as  follows.  (Note 
primes  denote  per  unit  length;  a = actual  value) 


Vbase 

Lbase 

C ' 

case 

Xbase 

Zbase 

^ase 

^base 

cbase 

Ibase 

tbase 

w, 

base 


VLN  (sVstem  peak) 


1^’  (L^  for  34>) 

C'  (Cl  for  34>) 

3 J. 

d (line  length) 


/L  ‘ /C ' 
a a 


Zbase/d 


L'd 

a 


Cad 


Vbase//Zbase 


✓ITcT  d 

a a 


i"J2qr a 


To  obtain  scaled  data,  divide  each  individual  para- 
meter by  its  base.  The  data  cards  are  as  follows: 
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1-  TMAX  , TPLOT  , DTPLOT  , LX  , NP  , IPLOT  (1<J)  & 3 <+> ) 

(F10.4)  (F10. 4 ) (F10.4)  (110)  (110)  (110) 


TMAX  - maximum  scaled  problem  time  for  program  to 
run;  at  least  2 cycles,  scaled,  for  a 
sinusoidal  source 

TPLOT  - scaled  problem  time  at  which  voltage  versus 
position  is  plotted.  Cannot  be  zero. 

DTPLOT  - scaled  problem  time  periods  after  TPLOT  at 
which  subsequent  plots  are  made.  Cannot  be 
zero. 

LX  - node  at  which  voltage  vs  time  is  plotted; 

1 < LX  < N=1 . 

NP  - determines  number  of  points  plotted  in  volt- 
age versus  time,  i.e.,  every  NP  points.  The 
calculating  At  is  fixed  internally  at  1/N 
(see  Section  2).  The  plotting,  and  printing, 
time  increment  is  NP*At. 

IPLOT  - plot  option 


1. 

voltage  versus 

position 

2. 

voltage  versus 

time 

3. 

both 

2- 

R 

L 

c 

N 

(140 

(F10.4) 

(F10.4) 

(F10 . 4 ) 

(HO) 

2- 

RA  , 

LA  , 

CA  , 

RB  , 

LB  , 

CB  , 

(340 

(F10.4) 

(F10. 4) 

( F10 . 4 ) 

(F10. 4) 

(F10.4) 

(F10.4) 

N 

(HO) 


R - scaled  line  resistance;  RA  and  RB  are  scaled 
zero  and  positive  sequence  values  respectively 

L - scaled  line  inductance;  LA  and  LB  are  scaled 
sequence  values 

C - scaled  line  capacitance;  CA  and  CB  are 
scaled  sequence  values 
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N - number  of  line  sections;  the  line  length 
divided  by  15  km  (9.3  mi)  to  the  nearest 
whole  number,  maximum  number  is  48 


EMAX  , OMEGA  , THETA  , RS  , LS 
(F10.4)  (F10.4)  (F10.4)  (F10.4)  (F10.4) 


EMAX  , OMEGA  , RSA  , LSA  , RSB  , RSC 
(F10.4)  (F10.4)  (F10.4)  (F10.4)  (F10.4)  (F10.4) 


EMAX  - maximum,  peak  value  (line  to  neutral  for 
3<j>)  of  the  voltage  source,  usually  1.0 


OMEGA  - 2 7T f 

THETA  - phase  shift  in  radians 

RS  - scaled  source  resistance;  RSA  and  RSB  are 
scaled  zero  and  positive  sequence  values 
respectively 

LS  - scaled  source  inductance;  LSA  and  LSA  are 
scaled  sequence  values 


TA  , TB  , TC  , THETA 
(F10.4)  (F10.4)  (F10.4)  (F10.4) 


TA  - scaled  time  delay  for  a-phase 

TB  - scaled  time  delay  for  b-phase 

TC  - scaled  time  delay  for  c-phase 

THETA  - phase  shift  in  radians 


GLL  , GAML  , CL 
(F10.4)  (F10.4)  (F10. 4) 


5-  GLLA  , GAML A , CLA  , GLLB  , GAMLB  , CLB 
(F10.4)  (F10.4)  (F10.4)  (F10.4)  (F10.4)  (F10.4) 


GLL  - scaled  load  conductance;  GLLA  and  GLLB  are 
scaled  zero  and  positive  sequence  values 
respectively 
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GAML  - scaled  load  gamma  (1/inductance) ; GAMLA  and 
GAMLB  are  scaled  sequence  values 

CL  - scaled  load  capacitance;  CLA  and  CLB  are 
scaled  sequence  values 

(Note;  the  load  can  be  any  parallel  combination  of 
inductance,  resistance,  and  capacitance.) 

If  a value  is  left  blank  on  a data  card  it  will  be  inter- 
preted as  zero  in  the  computer.  For  programs  that  run  for 
long  periods,  the  JCL  cards  controlling  run  time  may  have 
to  be  changed.  Always  check  the  last  data  cards  to  insure 
that  they  correspond  to  the  plot  options  chosen,  since  they 
label  the  plots. 
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Program  Listings 


/ J'JU  JOB , PAGE  5 = 60 « T 1 V |:  = 2 o o 

REAL  L,!S,lR,IL,llt,ICl,ILS,ICS,LS,IO 
D i M£fjS  I UN  XJ  ( 50)  . XL  AO  ( 5)  . V L A 11  ( 5 | , Gt  AR  ( r.  1 , 0 ATI  AB  ( 5 I 
DIMENSION  Y(50,50)  ,V  (501  , IB (50)  , PVI50)  ,PI  B(50  ).CM(50)  ,0(50,  501,  CC 
H SO)  ,CV(bO) , r ! ME  ( 900  1 , VOl.  T (900  I , AMPS  ( 900) 

RE AO (6 , 320  IT  HA X ■ TPIOT  .OTPLOT  , t X , NP.  i PL OT 

K t AH ( 5,3 3 0 ) ■ , L , C , N 

READ!  6, 360  )EM.\X, OMEGA, THETA, RS,tS  

RCA0(6.350)GU.,CL,GAML 

WRITC(6,360)  _ 

WRITE (6,370) 

WRITE  (6,38 OJ 

WRI  U(6,3  9 0_)  TMAX,  TPLOT  , OTPLOT,  LX,  NP,  I PLOT 

WRI  TE ( 6 , 600 ) _ 

WRI  TE  ( 6 , A 1 0 I 

WRITE(6,620)R,L,C,N 

WRITE (6,630) 

WRITE (6,660  ) 

WRIT E (6, 6 50) E MAX, OMEGA, THETA, RS,ES 

WR  ITE(6,660)  

WRITE (6,670) 

WRITEI6,680)GAML,CL,GLL  __  _ 

WR ITE(6,690) 

NC  = NP 

T = 0. 

DX=l./N 

N = N+  1 

DT = DX 

L=L*OX 

C=C*DX 

R=R*DX 

LS  = LS  + .000000l | 

NN=N+l 

_NPTS=(TMAX*(N-1  ) )/NP*l_  . 

DO  10  1=1 ,NN 

1 B C 1 )=0«  

V ( I ) = 0 • 

PV(I)=0.  

P1B  ( I 1=0. 

CM ( I 1=0. 

CC( I ) =0. 

CVII 1=0. 

DO  10  J= l , NN 

Y 1 1 , J ) =0 . 

10  B ( I , J ) =0. 

GG=1000000. 

GGS=l,/( (2.+ES/DT)*RS)  ' ' 

GL  = 0T*GAML/2.H  2. »CL/DT ) »GLl  

GS*l./( ( 2 . *L/DT  ) + R I 

GP=2.*C/DT  

C****BUILD  THE  Y-MATRIX**»* 

Y 1 1 » 1 ) = GG  *GGS  __  _ _ _ 

Y ( 1 , 2 1 =-GGS 

YJ?J  1 ) = - G G S 

DO  60  1=2, NN 

DO  50  J = 2 » NN  _ _ 

I F ( J . NE • I I GO  TO  30 

If II.WE.2IG0  TO  20  

Y ( I , J ) = Y ( I.JUGGS+GS 


CO  TO  *>0 

20  IFI l.NE.NNISO  TO  40  

Yll.JI=VII  , Jl+GS  + GP  + GL 

GO  TO  50  

30  1FIJ.N 6.1+1. AMO. J . NE • I - 1 ) GO  TO  50 

3ajaJJihl’.LL,.JJ.-iS 

CO  TO  50 

40  Y(  I , Jl-Ytl  ,J)+2.*GS  + GP  _ 

50  corn  i out 

C»**<MNVERT  Y-MATK IX*  **  + ..  

00  60  1=1, NN 

M-B-ll.  U-=g, 

00  70  1 = 1 ,11 

RAT  I 0=-Y  ( 1+1  , 1 )/Y(  I , I ) . 

YII  + 1 , 1 + 1 )=Y( I + l, I + 1 ) +RAT 1 0*Y ( 1 , 1+1 ) 

_ DO  70  J = 1 , ON  . _ 

IF(ABSIRATIO*B( I , Jl> .LT. 1E-101G0  TO  70 

B(  1 + 1 , -M=n  (1  + 1 , JURAT  I Q»B(  1 . J) 

70  CONTINUE 

DO  80  1 = 1,  N ..  _ 

K=NN-I 

. . DO  80  J=  l , NN  1 

RAT  I 0 = -Y ( K , K + ll/Y(K+l,K+ll 

I F ( A B S I R A T 1 0*  B ( K+l , J) I ,LT. IF-10)GO  TO  60 

B«K, J)  = B(K, J) +RAT 10*B(K  + 1 , J ) 

80  CONTINUE  ..  . 

DO  90  1=1 , NN 

DO  90  J = l , NN  _____  

90  6( I , J)=B(  I , JI/YI  I , I ) 

DO  1 I 0 I = 1 ,NN 

DO  110  J= 1 , NN 

J JF(ABS(B(I,J)).LT.1E-10JGO_TO  100_ 

YI I. J)  = B< I, J I 

CO  TO  110  __ 

100  Y 1 1 . J) =0.0 

110  CONTINUE 

VMAX--0. 

P I L 1=0. 

PICL=0. 

JX=1  ..  _ . . 

I F 1 1 PLOT . EQ . 0 ) GO  TO  130 

IF  I I PI  0T-2T1 ?0,  1 30, 1 70 

120  READ (5, 310)XLAB,YLAB,GLAB,DATLAB 

130  CONTINUE  _ 

IFIT.GT.THAXIGO  TO  260 

CC«1)=GGS*  IPV(  1 >-PV(2)l+PlBm*l  I2.*LS/DT  )-RS  ) / I l 2.*LS/DT  )+RS) 

DO  140  K=2,N 

140  CClKl =GS*I PV(K) -PVlK+1 ) ) + P I B ( K ) » ( ( 2 . *L /D T 1 — R ) / ( 1 2 . »l / DT ) +R ) 

DO  150  J = 3 « NN 

150  CVIJ)=PIB(  J-D-PIBI  J)+GP*PV(J)  

CLL=(DT*GAKL/2. 1 *PV 1 NN 1 +P  l U. 

CCL=-I2.*CL/0T ) *PV(NN)-P1CL  _ 

1 L=CL  L + CCL 

ANG=ORFGA*T+THETA 

1S=FKAX*C0SI ANGJ 

ig=  i snooooco.  

'CHI  1 )= IC-CCI 1 ) 

DO  160  J = 2,N  

160  CHI JI=CC( J-l J+CVt JI-CCI Jl 


80 


320  FORMA1 ( 3F 10.4, 3110) 

330  FORMA)  l 3T  10. '>.2110) 

340  FORMA  r ( sr  1 0.  4 ) 

350  F0f!MAT(3F10.4) 

360  TORMAT  (■!•,'*♦♦♦♦♦ 

1 1 l)f, »j* tmntn>Mi>ut<nm»n'  .//I 

370  FORMAT (44X, ' ♦♦♦♦OU1 PUT  CON TROL ♦*♦»',// ) 

380  FORMA  I ( 27X ,'THAX',6X,'TPL0T',4X,'DTPL0T'f8X,'LX',8X,,NP',6X,,lPL0T 

1 ■ ,//  1 

390  FORMAT (22X ,3F 1 0. 3 i 3 1 10 ) 

400  FORHAT I//.43X, • ’♦♦♦SCALED  LINE  OATA ♦♦♦♦',// ) 

_4JJ2_  _r  OR  MII25  A.'RISISTANCt'.ZX.  ' INDUCTANCE  ' ,2X.  ' C AP AC  1 TANCE  ' , 6X , 


I'SECT  IONS'  ,//> 

420  FORMAT (27X.F 10.3,2X,F10.3,3X,F10.3,2X,I10) 

430  FORMAT (// ,41 X ♦♦♦♦SC ALED  SOURCE  DATA****',//) 

44  0 FORMAT (31 X,' EM  AX' ,5X, 'OMEGA' ,5X . 'THETA' ,2X,  'RESISTANCE' . 
1 • INDUCTANCE ' , // 1 

450  FORMAT [ 2 5 X . 6 F 1 0 . 3 ) 


2X, 


460  FORMAT!//, 43X, •♦♦♦♦SCALED  LOAD  CA T A ♦♦♦♦',// 1 

470  FORMAT  I 39X , • GAMMA' , 2 X , • C A P AC  1 T ANC E ' ,2X,  'CONDUCTANCE ', // ) 

400  FORMAT (34X.F10. 3, 1X,F10.3,3X,F10. 3) 

490  FORMAT (///,*’ 

IT 

500  FORMAT ( 5X, ' VMAX^'  ,F 1 0. 3 , / / ) 

510  FORMAT (5X, 'NODE=* , 12,/ 1 

520  FORMAT  </,5X, « TIME* , 4 X , • VOLT  AGE • , 3X , ' CURRENT  •,  /) 

530  FORMAT (10F10.4) 

540  FORMAT (/,6X, 'NODE ' , 3X , • VOL T AGE • , 3X , 'CURRENT'  , / ) 

550  FORMAT!///, 5X, • TIME  *',F6.3) 

560  TORMAT  (IIP, SF  10. 4) . 


570 


FORMAT (/, • ' 

STOP 

END 


,/) 


'GO 


/DATA 


/DATA 
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/JOB  GROSS, PAGES=100,TIMC=300 

REAL) A, IB, IC. IS  A, ISO, ISC. IGA, IGB, !GC,LA,LB,LC, I RAX, I BA, IBB, IBC, 
1ILA, I LB , 1LC, 1LLA, ILLB, ILLC, 1CL A , ICLB, ICLC, II.SA, ILSB, ILSC, ICSA, 
2ICSB, 1CSC,LSA,LSB,LSC 

DIMENSION  YAI50.50) , Y8<50, 50), YC( 50,  50) , VA( 50) , V8( 50) , VCI50)  , 

l.l.!iA±5  I.L>  BULOLxJ  EC  I :■  OJ  J.SULA  (501  .PVBI50  U.R  Y C I MLUiU-E  A.L5.0)  ,£lg£.(5QJ-i 

2PI8CI 501 ,CMA( 50) ,CM3 ( 50) , CMC ( 50) , CC A ( 50 ) , CCB ( 50 ) , CCC ( 50 ) , CV A ( 50 ) , 

3CVB150) • C VC 1 5 0 ) .TIRE (500) , A VOL  I ( 500 > , BVOL T ( 500 ) , C VOl I 1 500), A AMPS 

5(500) , BAMPS (500  > .CAMPS! 500 ) , X J ( 50 ) , XL AB ( 5 ) , YL AB ( 5 ) ,GL AB ( 5 ) , 

5DAILABI5) .DATLAA15) .DATLACI5) 

READ! 5, 300) TMAX , TPLOT , UTPL 0T , L X , NP, I PLOT 

READ  1.5a_3  i 

READI5, 32 01EMAX, OMEGA, RSA,LSA,RSB,LSB 

. READ(5,320)TA,TB,TC, THETA  . 

REA0I5, 330  IGLL A, G AML  A, CL  A, GLLB, GAMLB.CLB 

. HRITE(6,350)  __1_  . . 

HR ITE(6,350) 

HR  I TC (6, 3 50) 

HR! TE (6,370 ) TMAX, TPLOT, DT PLOT, LX, NP, 1 PLOT 

HRITE  (6, 380)  : ..  _ 

HR) TE (6,390)N 

WRITE  (6,500)  .......  ...  . 

HRITE (6,510) 

HRITE(6,5?0)RA,LA.CA 

KRI TE ( 6 , 5 30  I 

HRITE(6,550)RB,LB,CB  . _ __  

HRI TE ( 6 , 550 ) 

HRI  TE  ( 6, 560)  __  . ... .....  

HRI TE (6,570  JEM AX, OMEGA, T A, TB.TC, THETA 

HR  1 T E ( 6 , 5 80  ) 

HRI TE (6,590) 

HR)  TE  (6,500)  RSA.LSA 

HRITEI6.510) 

HRITE(6,520)RSB,LSB  _ 

HRITE(6,530) 

HR) TE (6,550) 

HR) TE (6,550) 

HR!TE(6,560)GLLA,GAMLA,CLA  

HRITE (6,570) 

HRITE(6,580)GLLB, GAMLB.CLB  

HR  I TE ( 6 , 5 90 ) 

OX=)./H 

N = N*  1 

NN=N*  l _ _ 

NPTS-(1MAX*(N-1 ) l/NP+l 

DT=OX 

NC  = NP 

T = 0,  

(.  A = L A + DX 

LB  = LB*DX  

LC*IB 

CA=CA*OX  

CB=CB*OX 

CC^C 

RA=KA*OX 

_ RB  = RB*OX  _ 

RC  = RB 

RSC  = RSB  

LSA*LSA*. 0000001 
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istuistw  .OCOOOOl 

LSC=LSB 

GAMLC=GAMLB 

_ _ CLC^CLB  

GGA=iconooo. 

GC.R-GGA 

GGC=GGA 

GGSA-1./  ( [ 2. *LSA/UT  l+RSAI 

GLA  = ( DT*GAML A/2  . ) * ( 2.  *CL  A/DT  ) +GLLA 

GSA=l./{  (?.<‘LA/DT ) «RA)  . . ___  

GPA=2.*CA/DT 

C»»***»SUPROUTl g _TQ__!\U  1LD.  AND_-1N  VERT  THE  Y-.M  ATRIX 

CALL  YINVKT  < Y A,  GGA,GL  A,  GS  A,  CPA,  fiN,  GGSA  ) 

GGSli=l./l  (2.*LSB/0T)*RSB)  

GLB=  I DT*GAHLB/2.  ) + (2.0CLB/t)T)+3LLB 

GSB= 1 . / ( (2.fLB/DT)+RB)  

GPB=2.*CB/DT 

GGSC- GGSB 

GLC=GLB 

GSC  = GSB 

GPC=GPB 


CMB ( 1 1=0. 
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REA0(5,6A0I0ATLA8 

RE  AP  f 5 , GAOIDATLAC  _ 

30  CONTINUE 

1 F ( 1 .GT.THAX 1G0  TO  260 

CCA  ( I >=-GCSAMPVA(  1 ) -P  VA(  2 ) > +P  IB  At  1 ) * ( ( 2 . * L SA/O  T I -R  S A I / 

LLt2L.lLSm  J L+R3-AJ 

CCBtl)  =GG$B*  tPVBIH-PVBian  +PIBBI  n*(  (Z.*ISB/DT)-R$8)/ 
1 ( (2.*LSB/DT)*RSt>> 

CCC 1 1 )=GGSC*tPVC( 1 1-PVCI2  ) l+PIBCin  *( (2.*LSC/0T)-RSCI / 
1(  I2.*LSC/DT»*RSC) 

00  AO  K - 2 i N 

CCA(K)M.S.\*(PVAt<)-PVA(K+l)  I 4 P [ BA  ( K ) » ( ( 2 . »L  A/ DT) -R  A I / 

U (2.<  IA/0T)4RA) 

CCB(K)=GSB*(PVB(K)-PVBCK»1) J +P J BB t K ) ♦ t ( 2. *L6/DT ) -RB I / 

1 t 12.*LR/0T  )*RB> 

AO  CCCUMGSCMPVC  U)-PVC(K+l  I )*PIBC(K  )*t  (2.*LC/DT1-RC)/ 
l( (2.*LC/0T)tRC) 

00  50  JO.NN 

CVAI  J)=PI3A(  J-U-PIBAt  J I + GPA*PVA t J ) 

CVB(J)=PI88( J-l )-PlBB(J)+GPB*PVBU)  

50  CVCI J ) = P I BC I J-1)-PI8C( J ) 4GPC*PVC<  J ) 

CLLA=(DT*GAKLA/2.  ) *P  V A(  NN  ) + P I V.  I A 
CLL8=(0T*GAMLB/2. ) * P VB I NN I + P I LL B 

CH.C=tDl*GAMLC/?.  )»PVClNN)»PRtC 

CCLA=-(2.*CLA/DT)  *P V A ( NN ) -P 1 CL A 

CCtB=-(2.*CLB/0T)*PVB(NN)-PICLB 

CCLC=-(2.*CLC/DT)*PVC(NN)-PICLC 

ILA  = CLLA  + CCLA  

ILB=CLIB+CCLB 

ILCMLLCtCCLC 

1F(T-TA160,70,70 

60  I A = 0 . 1 

GO  TO  80 

70  ANGA  = OKfcGA*T  + THETA  

IA=EMAX*COS(  ANGA) 

BO  IFIT-TB)90.100,  100 

90  IB=0. 

GO  TO  110  

100  ANGB=0MEGA*T-2.09A 3951 ♦ THETA 

1 B=EMAX*COSl ANGB)  _ 

110  IF< T-TCI120,  130, 1 30 

120  1C-P. 

GO  TO  l AO 

_ 130  ANGC-OMCGA*T4?.09A395l+THFTA  _ 

IC=EMAX*COS< AN3C) 

•_1A0  CONTINUE  

ISAM  I At  I B+ IC ) / 3. 

I SB- ( 1 A- I B ) / 3 « 

ISCM  1A-  10/3. 

I GA  = I SAMOOOOOO.  _ _ 

IGB=ISB*1000000. 

IGC=ISC*1000000.  _ _ 

CMA 1 1 ) - I GA-CCA ( 1) 

CJHBJ  I ) = t Gfl-CCB  ( 1 ) 

CMC ( 1 1 * 1 6C-CCC ( 1 ) 

DO  150  J*2,N  . 

CMA(J)=CCA(J-1) +CVA ( J ) -CCA ( J ) 

CHB(JI-CCB(J-1 I *CVB(JI-CCB(Jl  

150  CMC ( J > =CCC  C J- 1 )*CVC(J)-CCC(J) 
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A 


CMA(NN)=CCA(NI *CVA(VN|-!LA 
CMB(NN)=CCBlN)*CVBl‘ri)-llB 
CMC(NN)=CCC( K) *CVC(NN)-!LC 
C **♦<"» tsuBRUUI  I\E  TO  CALCULATE  NUDE  VOLTAGES 
CALL  VXCMIVA,  YA.CMA.NR) 

CALL  VXC yjLT'gjJfii  ».CJ:L3  LU'll 

CALL  VXCMIVC.YC.CMC.NN) 

I OA ( 1 )=GGSAM  VA  (1 ) -V A ( 2 ) ) «CCA (1) 

I BB ( l )=GGSB*(VB( 1 )-VB(2) ) +CCB  ( 1 ) 

_IBC< l)=GGSC*I VC( 1)-VC(21) ♦CCC(l) 

DO  lf>0  K = 2,V 

LD A ( K 1 =GSA«(  VA(  K1-VA(  K<  1 ) 1 «CC M>;] 

1BB(K)  = GSB*( VB(K)-VB(K  + 1 ) I+CCBIK) 

160  IBC (KJ^GSCofVCTKl-VC  (K+l  ) 1 +CCC1K) 

IBA1 NN ) = GL A*V A ( UN  1 ♦ 1 L A 

IBBINN)=CLB*VB(NN)+1LB 

1BC(NN)=GLC*VCIN‘1)  + ILC 

1UA=(DT»GAMLA/?.  ) »VA  (N*n  «CLL  A 

1LLB={DT*GAMLB/2.)*VB(NM)*CLLB 

1LLC= (DTf&AKLC/2. ) *VC ( NR ) +CLLC 

ICLA-(2.*CLA/DT)*VA(NN)«CCLA 

ICLB= (2.0CLB/DT ) *VB (NN MCCLB 

ICLC=(2.*CLC/0T 1 *VC l NN 1 +CCLC 
C**»»*»»SIIBR0UT  U.'E  TO  CALCULATF  PHASE  VALUES 

CALL  M0T0PH(PVA,PVB,PVC,P1BA,P1BB.PIBC*NN) 

I F ( ABSI PVA (LXt l ) ) . LT. ABSI VAMAX) ) GO  TO  170 

VAMAX  = PVA<  LXt l ) 

1J0  IF(ABS(PV8(LX+1 ) I.LT.ABSIVBHAX) )GO  TO  180 

VBMAX^P VR ( LX* 1 ) 

180  1 F ( ABS I PV C I L X« 1 ) ) .LT.  ABS(  VCWAX1  1 GO  Tn  190 
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BVOLT(JX)=PVBlLX) 
CV0L1  (JX)=PVC(LX) 
AAMP$IJX)=P1«.'.(LX) 
B AMPS ( JX I -PICBtLX ) 
CAMPSIJX)-PIBC(LX) 

XJJ'.ULJXlil 

JX=JX*1 


2A0  NC=NC*l 

00  200  J=l,NN 
PVA ( J ) =VA ( J ) 
PVBl J)=VB( J) 

m-ijjjovx  uj 

P I B A ( J ) - I B A ( J ) 
PIBBI J)=IBB( J) 
250  P1BCI JJ-IBCI J) 
P1LLAMLLA 
PILLB=  ILLB 

■ . P1UC=IUC 

P 1 CL  A= I CL A 
P ICLB= ICLB 
P ICLC= ICLC 


T=r*DT 
GO  TO  30 
260  CONTINUE 


I F ( IPLOT.EO.OIGO  TO  290 
JF(IPLOT-2)290,270,270 
270  WRirE(6,6S0) 

WRI TE (6,690) 

WRITE(6,7001 VAMAX, VBMAX,  VCKAX 
WRITE(6,650)LX 


WR I TE ( 6 , 660 ) 

00  260  J = 1 ,NPTS 

28  0 HRITC(6,670|T!HE(I), AVOLT ( t ) ,BVOLT( I ) ,CVOLT (I ) , A AMPS t J ) , B AMPS  1 I ) ,C 

1 AMP S ( I ) 

READ! 5,6A0)XLA8, YLAB.GLAB.DATLAB 

CALL  GRAPH! NPTS .TIME, AVOLT, A, 107, 12. 0,8. 0,0. 0,0. 0,0. 0,-5.0,XLAB. 


I YLA8 , GL  A8 j OATL AB I 
READ(5,6A0)DATLAB 

CALL  GRAPHINPTS.T I ME , BVOL T , 9, 1 0 7 , 0. 0 , 8 . 0, 0.  0 , 0. 0 , 0 . 0 , -5 . 0 , XL AB , 

lYLABtGLABfDATLAB)  

REA0(5,6A0I0ATLAB 

CALL  GRAPH(NPTS,T|ME,CVOLT,8,107,0.0,8.0tO.O,0.0,0.0,-5.0iXLAB, 
lYLAB(GLABtOATLAB) 

290  CONTINUE  _ _ _ 

300  FORMAT (3FIO.A,3IIO) 

310  F0RMATI6F10.A.2I10)  _ 

320  FORMAT ( 8F 1 0. A ) 

330  FORMAT ( 6F 1 0. A ) 

350  FORMAT(AAX,'****OUTPUT  CONTROL****' ,//) 

360  FORMAT <27X, 'TMAX* ,6X, • TPIOT' ,AX, 'DTPLOT' ,8X, ' LX' ,8X, >NP' ,6X, • IPLOT 

l* ,// > 

370  FORMAT >?2X,3FIQ. 3. 31)0)  


380  FORMAT l // , A 3 X , ' ****SCALED  LINE  DATA****',//) 

390  FORMAT t A9Xt' SECT  I ONS=', 12,//) 

AOO  FORMAT ( 37 X, 'RES  I STANCE  * «2X, • I NDUC T ANCE • , 2 X, 'CAPACITANCE',//) 
A10  FORMAT 123X,'7ERO-SEQ' ) 

A20  FORMAT ( • ♦ • , 36X  , F l 0. A,2X,FI0.A,2X,F10.A,//) 


mu 


kV/r 


{ 


A06I624 
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680  FORMAT!//,  19X, 'VOLTAGE  MAXIMA',//) 

690_  FORMAT (15K, ' VA*  ,9X,«  VB'  ,9X, ' VC.*_»//1 

700  FORHATC9X,3F10.4> 

STOP T 

END 


» 


SUBROUTINE  Y I NVRT ( Y , GG, GL , GS, GP.NN, GGS > 

DIMENSION  Y(50,50),BI50,50) 

N=NN- l 

DO  10  1=1, NN 

PO-10  j ~ 1 1 N N 

YII,J)*0. 

10.  BII,J)=0. 

Yt 1 , 1 )=GG*GGS 

_Y(  l ,2  I =-GGS 

Y (2 , 1 ) =-GGS 

DO  SO  1=2, NN 

DO  50  J=2,NN 

IFIJ.NE.DGO  TO  30  

IFII.NE.2IG0  TO  20 

YtI,J)=YU,J)»GGS*GS 

GO  TO  50 


JLLLLiilL-i;  

Y 1 1 » J ) =Y ( I , J )*GS'GP+GL 


GO  TO  50  „ _ ' 

30  IFI  J.NE.Itl.ANO.J.NE.I-llGO  TO  50 


Y 1 1 , 0 ) *Y I 1 , J ) -GS 
GO  TO  50 


V 
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40  Ytl,J)=Y(l«J)«2.*GS+GP 

50  CONTINUE  

C**»*INVIRT  Y-MATR1X**** 

00  60  1 = 1.  NN  

60  B ( 1 , I ) = l . 

DQ-Jil-Ul^N 

RAT I0=-Y (1*1.1 I/YI  I , I ) 

Y ( I *1 , 1*1 1 = Yl I + 1 , Hl)tRATlQ*Y(!  ,I  + ll 

00  70  J=1 f NN 

IF(ABS(RATIO*BUf  jn.LT.lE-10)C0  TO  70 

B(I«1,J)=B(1*1,J>+RATI0*B(I,J> 

7.0-C0H.T.1NUE 

00  BO  1 = 1,  N 

K»NN-I  . . ._ _ ..  ... 

00  80  J= 1 , NN 

RAT10=-Y(K,K+l)/Y(K*l,K*ll 

IF(ABS(RAT10*B(K+1 , JIl.lT. IE-10IG0  TO  80 

B(K, J)=BIKtJ>»RAT10«BlK+l,J) 

80  CONTINUE 

DO  90  1 = 1, NN  . 

DO  90  J=l,NN 

90  B(I,J)=B(1,J)/YII.1> 

DO  110  I=l.NN 

DO  110  J=1 .NN 

1F< ABSIBI l , JII.LT. 1E-10JG0  TO  100 

YU,J)=BII,J)  . ....  _ 

GO  TO  110 

100  Y(  I , J I =0.0  _ •_ 

110  CONTINUE 

RETURN 

END 


SUBROUTINE  VXCM { v, y, CM, N I _ __ 

DIMENSION  VI501»Y(50,50)»CM150) 

DO  10  J=1,N 

V(J>=0. 

DO  10  K=l , N _ 

IFIABSIYI  J,KM>CM<K»  I .LT.IE-10IG0  TO  10 

V(Jl=V(J)*Y(J,K)*CM(K ) 

10  CONTINUE 

RETURN 

END 


SUBROUTINE  MOTOPHIVA.VB, VC, IBA, IBB, 18C,N) 

REAL  VA(50!,VB(50),VC(50) , IBAI50I , IBC ( 501 , I BBI50 I 
DIMENSION  PVA1 50) .PVBI50) ,PVC(50) ,PIBA(50) .PIBBI50) ,P1BC( 50) 
DO  10  J=1,N 

PVA(J)=VA(JI*VB(JI+VCIJ>  _ _ 

PIBA(JI=IBA(J)+IBBtJ)+IBC(J) 

PVBIJ)=VAI J>-2.*VB(JI+VC(J»  _ _ 

PIBBI J)=IBA( J)-2.*IBB( JI+IBCI Jl 

PVCI JI=VA( JI«VB( J)-2.»VC(  J) 

10  Pl8:i  JI  = I6A(  J|*IBBUI-2.*IBCIJI 

_ 00  20  J=1,N  _ 

VA(JI=PVA(J>  * ~ ----- 

VB( JI=PV8( Jl 

vcui=pvcij>  ""  “ ‘ ~ 


